#!/usr/bin/perl
use Getopt::Std;
#getopt("a:c:d:n:m", \%args);
#getopt("a:d:n:m", \%args);
$ambi = "chr8q24_GBM.tped";
chomp $ambi;
#$chr = $args{c};
#chomp $chr;
$dosage = "impute_8q24_out_ambi_final";
chomp $dosage;
$newdosage = "chr8q24_GBM_2step_ambi.tped";
chomp $newdosage;
if($ambi eq "")
{
        die "entered ambi snps file name  is empty\n";
}
#if($chr eq "")
#{
#        die "entered  $chr  is empty\n";
#}
if($dosage eq "")
{
        die "entered dosage file name is empty\n";
}
#$grep = "grep -P '^".$chr." ' $ambi |";
if(!(-e $ambi))
{
	die "$ambi file not exists\n";
}
#$grep = "gunzip -c $ambi|grep -P '^".$chr." ' |";
#die $grep."\n";
#open(SAM,"$grep") or die "no file found $ambi\n";
#open(SAM,"gunzip -c $ambi|") or die "no file found $ambi\n";
#open SAM,"grep -P $ambi" ;
open(SAM,"$ambi") or die "no file found $ambi\n";
open(DOS,"$dosage") or die "no file found $dosage\n";
open(WRBUFF1,">$newdosage") or die "unable to write to ambi dosage file\n"; 
#$line = <DOS>;

#reading all the ambiguous snp positions in to hash

while($line = <SAM>)
{
	chomp $line;
	$line =~ s/\n//g;
	$line =~ s/\r//g;
	@array_sam =split(/\s/,$line);
	if(($line =~ m/ T/ && $line =~ m/ A/) ||($line =~ m/ G/ && $line =~ m/ C/))
	{
		$hash_sam{$array_sam[3]} =$line;
		#print $line."\n";
	}
}
undef(@array_sam);
$num_samp = 0;
#reading the dosage file
while($line = <DOS>)
{	
	chomp($line);
	@array = split(" ",$line);
		if(exists($hash_sam{$array[0]}))
		{
			#print "$line\n";
			$line1 = $hash_sam{$array[0]};
			@array_sam = split(/\s/,$line1);
			$a1 = $array[2];
			$a2 = $array[3];
			$count_fwd = 0;
			$count_rev = 0;
			shift(@array_sam);
			shift(@array_sam);
			shift(@array_sam);
			$pp1=shift(@array_sam);
			$line1 = join(" ",@array_sam);
			$j = 0;
			#print "$line\n$line1\n";
			#count the maf and checking should be greater than 30%
			$maf_check = $line1;
			
			
				#print "counting the matches of homozygotes in both strands\n";
				for($i=4;$i<@array;$i++)
				{
					$d1 = $array[$i++];
					$d2 = $array[$i++];
					$d3 = $array[$i];
					$gnt1 = $array_sam[$j++]." ".$array_sam[$j++];
					$gnt1_flip = $gnt1;
					$gnt1_flip =~ tr/ATGC/TACG/;
					if($d1 >0.5)
					{
						$gnt2 = "$a1 $a1";
					}
					elsif($d2 > 0.5)
					{
						$gnt2 ="$a1 $a2";
					}
					elsif($d3 > 0.5)
					{
						$gnt2 = "$a2 $a2";
					}
					else
					{
						$gnt2 = "0 0";
					}
					@g1 = split(" ",$gnt1);
					@g1_f = split(" ",$gnt1_flip);
					@g2 = split(" ",$gnt2);
					#die  "$pp1 gnt1:$gnt1\t@g1\t\t\t\tgnt-flip:$gnt1_flip\t@g1_f\t\t\t\t\tgnt2:$gnt2\t@g2\n";
					#print "$gnt1\n";
					if($g1[0] eq $g1[1] && $gnt1 ne "0 0" && $g2[0] eq $g2[1] && $gnt2 ne "0 0" && $gnt1 eq $gnt2)
					{
						$count_fwd++;
					}
					if($g1_f[0] eq $g1_f[1] && $gnt1_flip ne "0 0" && $g2[0] eq $g2[1] && $gnt2 ne "0 0" && $gnt1_flip eq $gnt2)
					{
						$count_rev++;
					}
				}
				undef(@array_sam);
				#print "$count_fwd $count_rev $array[0] $array[1] $array[2] $array[3] $array[4]\n";
				#fliping the input genotype according to counts
				if($count_fwd > $count_rev)
				{
					#$finalgnt = $line1;
					#die  "$count_fwd $count_rev $array[0] $array[1] $array[2] $array[3] $array[4]\n";
				}
				else
				{
					$flip{$array[0]}=1
					#$line1 =~ tr/ATGC/TACG/;
					#$finalgnt = $line1;
					#die  "$array[0] $array[1] $array[2] $array[3] $array[4]\n";
				}
		}		
}
open(SAM,"$ambi") or die "no file found $ambi\n";
while(<SAM>)
{
	chomp($_);
	@array=split(" ",$_);
	if(exists($flip{$array[3]}))
	{
		print WRBUFF1 shift(@array);
		print WRBUFF1 " ".shift(@array);
		print WRBUFF1 " ".shift(@array);
		print WRBUFF1 " ".shift(@array);
		$_ =join(" ",@array);
		$_ =~ tr/ATGC/TACG/;
		print WRBUFF1 " $_\n";	
	}
	else
	{
		print WRBUFF1 $_."\n";
	}		
}
